US COVID Data

Data Prep

us <- read.csv("https://raw.githubusercontent.com/JaclynCoate/6373_Time_Series/master/TermProject/Data/USdaily.7.26.20.csv", header = T, strip.white = T)
us <- transform(us, date = as.Date(as.character(date), "%Y%m%d"))
us <- subset(us, select = -c(states, dateChecked, hospitalized, lastModified, total, posNeg, totalTestResultsIncrease, hash))
us[is.na(us)] <- 0
#Selecting only those dates with reported current hospitilizations
#us <- dplyr::slice(us,1:132)
us = us[order(as.Date(us$date, format = "%Y%m%d")),]
head(us)
##           date positive negative pending hospitalizedCurrently
## 187 2020-01-22        2        0       0                     0
## 186 2020-01-23        2        0       0                     0
## 185 2020-01-24        2        0       0                     0
## 184 2020-01-25        2        0       0                     0
## 183 2020-01-26        2        0       0                     0
## 182 2020-01-27        2        0       0                     0
##     hospitalizedCumulative inIcuCurrently inIcuCumulative
## 187                      0              0               0
## 186                      0              0               0
## 185                      0              0               0
## 184                      0              0               0
## 183                      0              0               0
## 182                      0              0               0
##     onVentilatorCurrently onVentilatorCumulative recovered death
## 187                     0                      0         0     0
## 186                     0                      0         0     0
## 185                     0                      0         0     0
## 184                     0                      0         0     0
## 183                     0                      0         0     0
## 182                     0                      0         0     0
##     totalTestResults deathIncrease hospitalizedIncrease negativeIncrease
## 187                2             0                    0                0
## 186                2             0                    0                0
## 185                2             0                    0                0
## 184                2             0                    0                0
## 183                2             0                    0                0
## 182                2             0                    0                0
##     positiveIncrease
## 187                0
## 186                0
## 185                0
## 184                0
## 183                0
## 182                0

Current Hospitalization Realization

Traits:

  • Heavy wandering behavior
  • What appears to be some noise that could be pseudo-cyclic behavior hidden by the large numbers.
ggplot(data = us, aes(x=date, y=hospitalizedCurrently))+
  geom_line(color="orange")+
  labs(title = "Current COVID Hospitalized Cases US", y = "Thousands", x = "") +
    theme_fivethirtyeight()

Stationarity

It is difficult to assume stationarity for this data due to multiple factors. We are working under the assumption that COVID is a novel virus and cases as well as hospitalizations will eventually return to zero. This being said our current modeling techniques do things such as return to the median or mimic the previously seen trends. Also, we see a heavy wandering trend in both new cases and hospitalization would be dependent on this as well as time. We will review the data and see what, if any, non-stationary components reveal themselves and model the data accordingly.

Forecast Current Hospitalizations

It is difficult to assume stationarity for this data due to multiple factors. We are working under the assumption that COVID is a novel virus and cases as well as hospitalizations will eventually return to zero. This being said our current modeling techniques do things such as return to the median or mimic the previously seen trends. Also, we see a heavy wandering trend in both new cases and hospitalization would be dependent on this as well as time. We will review the data and see what, if any, non-stationary components reveal themselves and model the data accordingly.

Variable Review

When reviewing the variables and the scatter plot from our previous EDA we can see that there are some correlations that were expected, for example currentHospitalizations is correlated with the variables that refelct ICU and Ventilaor patients. These metrics wouldn’t exist without hospitalizations. These variables also cannot help up predict hospitalizations because they after occurances. So we will actually be leveraing variables such as positive and positive increase in order to see if there is some sort of correlation between hospitalized patients and the number of positive cases.

  • Positive Tests Trend
ggplot(data = us, aes(x=date, y=positive))+
  geom_line(color="orange")+
  labs(title = "Total Cases US", y = "Millions", x = "") +
    theme_fivethirtyeight()

  • Positive Increase Trend
ggplot(data = us, aes(x=date, y=positiveIncrease))+
  geom_line(color="orange")+
  labs(title = "Increase in COVID Cases US", y = "Thousands", x = "") +
    theme_fivethirtyeight()

MLR w/ Correlated Errors for Currently Hospitalized Patients

Forecast Independent Variables

  1. Forecast COVID Cases
    1. Evaluation: Heavy wandering, frequency peak at 0, and slowly dampening ACF consistend with (1-B) component
#a
plotts.sample.wge(us$positive)

## $autplt
##  [1] 1.0000000 0.9802449 0.9605622 0.9410048 0.9216948 0.9025967 0.8836977
##  [8] 0.8649429 0.8462650 0.8277641 0.8094575 0.7914950 0.7738194 0.7563783
## [15] 0.7391523 0.7221011 0.7052654 0.6886822 0.6724063 0.6563621 0.6405962
## [22] 0.6250004 0.6094678 0.5940546 0.5788727 0.5639282
## 
## $freq
##  [1] 0.005347594 0.010695187 0.016042781 0.021390374 0.026737968
##  [6] 0.032085561 0.037433155 0.042780749 0.048128342 0.053475936
## [11] 0.058823529 0.064171123 0.069518717 0.074866310 0.080213904
## [16] 0.085561497 0.090909091 0.096256684 0.101604278 0.106951872
## [21] 0.112299465 0.117647059 0.122994652 0.128342246 0.133689840
## [26] 0.139037433 0.144385027 0.149732620 0.155080214 0.160427807
## [31] 0.165775401 0.171122995 0.176470588 0.181818182 0.187165775
## [36] 0.192513369 0.197860963 0.203208556 0.208556150 0.213903743
## [41] 0.219251337 0.224598930 0.229946524 0.235294118 0.240641711
## [46] 0.245989305 0.251336898 0.256684492 0.262032086 0.267379679
## [51] 0.272727273 0.278074866 0.283422460 0.288770053 0.294117647
## [56] 0.299465241 0.304812834 0.310160428 0.315508021 0.320855615
## [61] 0.326203209 0.331550802 0.336898396 0.342245989 0.347593583
## [66] 0.352941176 0.358288770 0.363636364 0.368983957 0.374331551
## [71] 0.379679144 0.385026738 0.390374332 0.395721925 0.401069519
## [76] 0.406417112 0.411764706 0.417112299 0.422459893 0.427807487
## [81] 0.433155080 0.438502674 0.443850267 0.449197861 0.454545455
## [86] 0.459893048 0.465240642 0.470588235 0.475935829 0.481283422
## [91] 0.486631016 0.491978610 0.497326203
## 
## $db
##  [1]  17.7475285  10.1272121   8.1941005   6.1372622   3.6673423
##  [6]   2.1381890   0.8844997  -0.3857354  -1.2308930  -2.2716130
## [11]  -3.0764132  -3.8707782  -4.5992105  -5.2234968  -5.7429794
## [16]  -6.3209084  -6.7999549  -7.2711475  -7.7287507  -8.2150804
## [21]  -8.5822471  -9.0368262  -9.3274696  -9.7150034  -9.8793477
## [26] -10.3026694 -10.2171729 -10.9486831 -11.1792449 -11.5414722
## [31] -11.7533685 -12.0032382 -12.2398160 -12.4967820 -12.6866927
## [36] -12.9359621 -13.1356344 -13.3714205 -13.5407817 -13.7677336
## [41] -13.9011781 -14.0623013 -14.1947429 -14.4218132 -14.5035227
## [46] -14.7049957 -14.8554851 -14.9990579 -15.1602596 -15.3080424
## [51] -15.3983431 -15.5036313 -15.5477357 -15.7601162 -15.8989007
## [56] -16.0472662 -16.1139387 -16.1931787 -16.3168599 -16.3675392
## [61] -16.4321318 -16.5678159 -16.5664973 -16.7389082 -16.8089441
## [66] -16.8535974 -16.9862537 -17.0350263 -17.1077034 -17.1394212
## [71] -17.2054793 -17.2719415 -17.3213157 -17.3792346 -17.4074767
## [76] -17.4795400 -17.5155214 -17.5405625 -17.5571892 -17.5626813
## [81] -17.6100773 -17.6142828 -17.6551199 -17.6909035 -17.7031509
## [86] -17.7473773 -17.7644017 -17.7877767 -17.8298413 -17.8244288
## [91] -17.8613989 -17.8180327 -17.7990496
## 
## $dbz
##  [1]  12.3783835  11.9824059  11.3184607  10.3811694   9.1646548
##  [6]   7.6655606   5.8897118   3.8659261   1.6709286  -0.5403053
## [11]  -2.5411473  -4.1440453  -5.3411027  -6.2600784  -6.9959398
## [16]  -7.5519535  -7.9149366  -8.1399736  -8.3445616  -8.6412919
## [21]  -9.0890880  -9.6819745 -10.3587036 -11.0280982 -11.6114762
## [26] -12.0816948 -12.4632524 -12.7940444 -13.0899873 -13.3416750
## [31] -13.5379862 -13.6909936 -13.8385604 -14.0243653 -14.2742680
## [36] -14.5843516 -14.9242721 -15.2532895 -15.5412795 -15.7814483
## [41] -15.9855807 -16.1678531 -16.3324279 -16.4740452 -16.5887688
## [46] -16.6841796 -16.7796466 -16.8965974 -17.0467154 -17.2259573
## [51] -17.4173594 -17.6006304 -17.7626750 -17.9021554 -18.0255117
## [56] -18.1387710 -18.2422649 -18.3320077 -18.4057674 -18.4682177
## [61] -18.5305656 -18.6048207 -18.6970293 -18.8040491 -18.9156848
## [66] -19.0206480 -19.1124394 -19.1913149 -19.2614785 -19.3263109
## [71] -19.3855708 -19.4364663 -19.4772892 -19.5103587 -19.5417050
## [76] -19.5776554 -19.6209105 -19.6689409 -19.7157844 -19.7559847
## [81] -19.7879270 -19.8142648 -19.8393453 -19.8657639 -19.8926241
## [86] -19.9166161 -19.9349149 -19.9475839 -19.9575982 -19.9684882
## [91] -19.9814602 -19.9941979 -20.0022876
b. Differencing Data: still observe a slowly dampening ACF, secondary (1-B) transformation
#b
pos.diff = artrans.wge(us$positive, phi.tr = 1, lag.max = 100)

c. Differencing Data: much more stationary data and have surfaced a seasonaly component = 7 seen on ACF peaks at 7, 14, 21

#c
pos.diff.2 = artrans.wge(pos.diff, phi.tr = 1, lag.max = 100)

d. Seaonality Transformation: stationary data, and an ACF that reflects data other than whitenoise to be modeled.

#d
pos.diff.seas = artrans.wge(pos.diff.2,phi.tr = c(0,0,0,0,0,0,1), lag.max = 100)

e. Model IDing: diagnose with aic.wge to determine phi, thetas: Both AIC/BIC select ARMA(4,2)

#e
aic5.wge(pos.diff.seas)
## ---------WORKING... PLEASE WAIT... 
## 
## 
## Five Smallest Values of  aic
##       p    q        aic
## 15    4    2   15.62972
## 9     2    2   15.75292
## 18    5    2   15.75381
## 12    3    2   15.75949
## 7     2    0   15.76130
aic5.wge(pos.diff.seas, type = "bic")
## ---------WORKING... PLEASE WAIT... 
## 
## 
## Five Smallest Values of  bic
##       p    q        bic
## 15    4    2   15.75485
## 2     0    1   15.79787
## 7     2    0   15.81493
## 4     1    0   15.82044
## 3     0    2   15.82122
f. White Noise Test: Reject the H_null and accept this data is not white noise
#f
acf(pos.diff.seas)

ljung.wge(pos.diff.seas)$pval
## Obs -0.3711871 -0.02127298 0.07854863 0.05194626 -0.1044837 0.3276797 -0.4658846 0.1565308 0.1284193 -0.1060784 -0.07693315 0.103073 -0.06391544 0.05002948 0.004471733 -0.03640371 -0.01551969 0.04133791 -0.0111003 -0.08153887 0.06570589 -0.03709668 -0.02390002 0.04588505
## [1] 1.568856e-12
#ljung.wge(pos.diff.seas, K = 48)$pval
g. Estiamte Phis, Thetas: 
#g
est.pos.seas = est.arma.wge(pos.diff.seas, p = 4, q = 2)
## 
## Coefficients of Original polynomial:  
## -1.8726 -1.4973 -0.4734 0.0283 
## 
## Factor                 Roots                Abs Recip    System Freq 
## 1+1.1292B+0.6983B^2   -0.8086+-0.8822i      0.8356       0.3681
## 1+0.7944B             -1.2588               0.7944       0.5000
## 1-0.0510B              19.6081               0.0510       0.0000
##   
## 
mean(us$positive)
## [1] 1175047
h. Model Building
i. Forecast 7/90 Days 
#h
#Equation building
#i
pos.preds7 = fore.aruma.wge(us$positive, phi = est.pos.seas$phi, theta = est.pos.seas$theta, d=2, s=7, n.ahead = 7)

pos.preds90 = fore.aruma.wge(us$positive, phi = est.pos.seas$phi, theta = est.pos.seas$theta, d=2, s=7, n.ahead = 90)

i. Plotting forecasts

plot(seq(1,187,1), us$positive, type = "l", xlim = c(0,195), ylim = c(0,5000000), ylab = "Positive COVID Cases", main = "7 Day COVID Cases Forecast")
lines(seq(188, 194,1), pos.preds7$f, type = "l", col = "red")

plot(seq(1,187,1), us$positive, type = "l", xlim = c(0,280), ylim = c(0,9000000), ylab = "Positive COVID Cases", main = "90 Day COVID Cases Forecast")
lines(seq(188, 277,1), pos.preds90$f, type = "l", col = "red")

  1. Forecast Increase Positive Cases
    1. Evaluation: Slowly dampening ACF and heavy wandering consistent with a (1-B) factor. As well as frequency peaks at 0 and .14. Consistent with (1-B) and seasonailty component of 7.
#a
plotts.sample.wge(us$positiveIncrease)

## $autplt
##  [1] 1.0000000 0.9740745 0.9439116 0.9109437 0.8869924 0.8711719 0.8605261
##  [8] 0.8462935 0.8144161 0.7772852 0.7356766 0.7045767 0.6829609 0.6646632
## [15] 0.6432819 0.6088202 0.5688872 0.5293417 0.5001353 0.4761589 0.4592724
## [22] 0.4426271 0.4148783 0.3789333 0.3458221 0.3182895
## 
## $freq
##  [1] 0.005347594 0.010695187 0.016042781 0.021390374 0.026737968
##  [6] 0.032085561 0.037433155 0.042780749 0.048128342 0.053475936
## [11] 0.058823529 0.064171123 0.069518717 0.074866310 0.080213904
## [16] 0.085561497 0.090909091 0.096256684 0.101604278 0.106951872
## [21] 0.112299465 0.117647059 0.122994652 0.128342246 0.133689840
## [26] 0.139037433 0.144385027 0.149732620 0.155080214 0.160427807
## [31] 0.165775401 0.171122995 0.176470588 0.181818182 0.187165775
## [36] 0.192513369 0.197860963 0.203208556 0.208556150 0.213903743
## [41] 0.219251337 0.224598930 0.229946524 0.235294118 0.240641711
## [46] 0.245989305 0.251336898 0.256684492 0.262032086 0.267379679
## [51] 0.272727273 0.278074866 0.283422460 0.288770053 0.294117647
## [56] 0.299465241 0.304812834 0.310160428 0.315508021 0.320855615
## [61] 0.326203209 0.331550802 0.336898396 0.342245989 0.347593583
## [66] 0.352941176 0.358288770 0.363636364 0.368983957 0.374331551
## [71] 0.379679144 0.385026738 0.390374332 0.395721925 0.401069519
## [76] 0.406417112 0.411764706 0.417112299 0.422459893 0.427807487
## [81] 0.433155080 0.438502674 0.443850267 0.449197861 0.454545455
## [86] 0.459893048 0.465240642 0.470588235 0.475935829 0.481283422
## [91] 0.486631016 0.491978610 0.497326203
## 
## $db
##  [1]  14.3997611  15.8955909   7.7642205   8.7050186   3.8198773
##  [6]  -1.1863576   2.6815084  -0.4667544  -1.0582274  -3.9729372
## [11]  -4.4388886  -5.3739230  -7.0574666  -3.9001118  -5.3732547
## [16]  -5.9203161  -6.0265835  -6.3073466  -8.7427682  -8.5315257
## [21]  -8.1694048  -8.8726978  -6.2228751  -5.6255282  -5.7869081
## [26]  -2.3190255  -0.6895940 -18.2510780 -12.7252753 -17.9895004
## [31] -13.1495229 -17.1809062 -12.9879824 -18.5885318 -14.0950716
## [36] -19.2823521 -16.1802599 -22.2193540 -17.1778910 -15.9864725
## [41] -13.8720977 -14.3606703 -12.3482825 -19.0159734 -12.8299009
## [46] -16.1634117 -19.2402721 -16.9940391 -19.9193164 -24.4462952
## [51] -12.2184745 -14.1347405 -11.1407485 -19.5608640 -22.3406178
## [56] -26.2948840 -15.9219942 -19.7896219 -16.3429699 -16.6545137
## [61] -13.2991307 -20.9465604 -12.5392646 -23.3015328 -21.5282804
## [66] -17.4748203 -26.6390880 -24.5923917 -17.5647166 -19.7282494
## [71] -17.6845420 -23.0825883 -19.0743556 -19.4958696 -19.5138333
## [76] -20.4270518 -18.7343253 -16.8982195 -14.1306252 -14.2054286
## [81] -15.0466299 -13.4474119 -14.7662082 -15.2848488 -14.7767269
## [86] -17.2543647 -17.5524024 -16.7680943 -23.6931393 -19.6888637
## [91] -22.9328710 -17.1805763 -15.7422487
## 
## $dbz
##  [1]  12.1744244  11.8068871  11.1912920  10.3234897   9.1987290
##  [6]   7.8134517   6.1690554   4.2794940   2.1855400  -0.0227779
## [11]  -2.1856123  -4.0930116  -5.5814407  -6.6274820  -7.3078903
## [16]  -7.6988921  -7.8526856  -7.8271468  -7.6919809  -7.5061868
## [21]  -7.3030888  -7.0977786  -6.9055053  -6.7545171  -6.6859448
## [26]  -6.7448940  -6.9710578  -7.3933729  -8.0283548  -8.8797888
## [31]  -9.9376266 -11.1745795 -12.5395304 -13.9486256 -15.2804168
## [36] -16.3913300 -17.1665125 -17.5812779 -17.7097183 -17.6695854
## [41] -17.5637030 -17.4562894 -17.3753921 -17.3239844 -17.2913960
## [46] -17.2628022 -17.2262162 -17.1767997 -17.1186074 -17.0640034
## [51] -17.0309217 -17.0382729 -17.1002705 -17.2209453 -17.3905128
## [56] -17.5857017 -17.7761429 -17.9367144 -18.0606665 -18.1649087
## [61] -18.2832201 -18.4521578 -18.6979308 -19.0286543 -19.4319317
## [66] -19.8765383 -20.3176913 -20.7053075 -20.9929698 -21.1442791
## [71] -21.1356798 -20.9585908 -20.6228760 -20.1586258 -19.6118560
## [76] -19.0349178 -18.4769322 -17.9782467 -17.5690679 -17.2703544
## [81] -17.0951532 -17.0493949 -17.1318377 -17.3332302 -17.6350503
## [86] -18.0085946 -18.4157858 -18.8133902 -19.1613062 -19.4324915
## [91] -19.6189508 -19.7292633 -19.7788462
b. Differencing Data: much more stationary data and have surfaced a seasonaly component = 7 seen on ACF peaks at 7, 14, 21
#b
inpos.diff = artrans.wge(us$positiveIncrease, phi.tr = 1, lag.max = 100)

c. Seaonality Transformation: stationary data, and an ACF that reflects data other than whitenoise to be modeled.

#c
inpos.diff.seas = artrans.wge(inpos.diff,phi.tr = c(0,0,0,0,0,0,1), lag.max = 100)

d. Model IDing: diagnose with aic.wge to determine phi, thetas: Both AIC/BIC select ARMA(4,2)

#d
aic5.wge(inpos.diff.seas)
## ---------WORKING... PLEASE WAIT... 
## 
## 
## Five Smallest Values of  aic
##       p    q        aic
## 15    4    2   15.62368
## 9     2    2   15.74701
## 18    5    2   15.74771
## 12    3    2   15.75351
## 7     2    0   15.75551
aic5.wge(inpos.diff.seas, type = "bic")
## ---------WORKING... PLEASE WAIT... 
## 
## 
## Five Smallest Values of  bic
##       p    q        bic
## 15    4    2   15.74832
## 2     0    1   15.79201
## 7     2    0   15.80893
## 4     1    0   15.81457
## 3     0    2   15.81522
e. White Noise Test: Reject the H_null and accept this data is not white noise
#e
acf(inpos.diff.seas)

ljung.wge(inpos.diff.seas)$pval
## Obs -0.3711871 -0.02127283 0.07854862 0.05194601 -0.1044838 0.3276797 -0.4658847 0.1565306 0.1284186 -0.1060791 -0.07693327 0.1030723 -0.06391648 0.05002847 0.004471288 -0.03640442 -0.01551996 0.04133725 -0.01110058 -0.08153895 0.06570572 -0.03709739 -0.02390062 0.04588417
## [1] 1.246447e-12
f. Estiamte Phis, Thetas: 
g. Model Building
#f
est.inpos.seas = est.arma.wge(inpos.diff.seas, p = 4, q = 2)
## 
## Coefficients of Original polynomial:  
## -1.8727 -1.4974 -0.4735 0.0283 
## 
## Factor                 Roots                Abs Recip    System Freq 
## 1+1.1293B+0.6983B^2   -0.8086+-0.8822i      0.8357       0.3681
## 1+0.7944B             -1.2588               0.7944       0.5000
## 1-0.0510B              19.6244               0.0510       0.0000
##   
## 
mean(us$positiveIncrease)
## [1] 22567.12
#g
h. Forecast 7/90 Days 
#7 day
inpos.preds7 = fore.aruma.wge(us$positiveIncrease, phi = est.inpos.seas$phi, theta = est.inpos.seas$theta, d=1, s=7, n.ahead = 7)

#90 day
inpos.preds90 = fore.aruma.wge(us$positiveIncrease, phi = est.inpos.seas$phi, theta = est.inpos.seas$theta, d=1, s=7, n.ahead = 90)

i. Plotting forecasts

#7 day forecast
plot(seq(1,187,1), us$positiveIncrease, type = "l", xlim = c(0,195), ylim = c(0,80000), ylab = "Increase in COVID Cases", main = "7 Day Increase in COVID Cases Forecast")
lines(seq(188, 194,1), inpos.preds7$f, type = "l", col = "red")

#90 day forecast
plot(seq(1,187,1), us$positiveIncrease, type = "l", xlim = c(0,280), ylim = c(0,80000), ylab = "Increase in COVID Cases", main = "90 Day Increase in COVID Cases Forecast")
lines(seq(188, 277,1), inpos.preds90$f, type = "l", col = "red")

Forecast Currently Hospitalized Patients using Positive Increase and Total Cases

MLR w/ Correlated Errors w/ No Lagged Variables
  1. Fit simple regression model predicting hospitalizedCurrently with positive, positiveIncrease and a trend.
us.mlr.fit = lm(hospitalizedCurrently~positiveIncrease+positive+date, data=us)
summary(us.mlr.fit)
## 
## Call:
## lm(formula = hospitalizedCurrently ~ positiveIncrease + positive + 
##     date, data = us)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -21825.6  -9253.8   -426.3   8552.6  24559.2 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -6.090e+06  9.915e+05  -6.142 4.96e-09 ***
## positiveIncrease  9.013e-01  9.757e-02   9.237  < 2e-16 ***
## positive         -1.452e-02  2.546e-03  -5.703 4.64e-08 ***
## date              3.327e+02  5.411e+01   6.149 4.76e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 12220 on 183 degrees of freedom
## Multiple R-squared:  0.7136, Adjusted R-squared:  0.7089 
## F-statistic:   152 on 3 and 183 DF,  p-value: < 2.2e-16
  1. Diagnose with aic.wge and see what the p, q would be with residuals from above model
    1. Below we can see that our aic.wge function has selected an AR(6) model for modeling our cmort information.
us.mlr.phi= aic.wge(us.mlr.fit$residuals, p = 0:12, q = 0:12)
us.mlr.phi
## $type
## [1] "aic"
## 
## $value
## [1] 15.33178
## 
## $p
## [1] 8
## 
## $q
## [1] 6
## 
## $phi
## [1]  1.42405076 -0.24799231 -0.61685924  0.06377608  0.97333749 -0.77347228
## [7]  0.41766520 -0.27190429
## 
## $theta
## [1]  0.6015312  0.1219283 -0.4916272 -0.4838232  0.7755189 -0.4855814
## 
## $vara
## [1] 3880013
  1. Now forcast with arima function with phi’s from above coefficients and AR(6) model #Error, can’t seem to figure out what’s going on
#us.mlr.modelfit = arima(us$hospitalizedCurrently, order = c(8,0,6), xreg = cbind(us$positive, us$hospitalizedIncrease, us$date))
#AIC(us.mlr.modelfit)

#fit = arima(CM$cmort,order = c(phi$p,0,phi$q), seasonal = list(order = c(1,0,0), period = 52), xreg = cbind(CM$temp, CM$part, CM$Week))

#AIC(fit) #3166
  1. Now we check the residuals for white noise from our above model
#acf(us.mlr.modelfit$residuals) #Visual ACF test
#us.mlr.wn.test = ljung.wge(us.mlr.modelfit$residuals)
#us.mlr.wn.test$pval

With the inability to reject the N null for the white noise test for an MLR model w/ correlated errors we will have to settle on the model being white noise with a seasonal and (1-B) component.

  1. Load the forecasted positive and increase positive in a data frame
#7 forecast
#regressions7 = data.frame(positiveIncrease = inpos.preds7$f, positive = pos.preds7$f, date = seq(188,194,1))
#90 forecast
#regressions90 = data.frame(positiveIncrease = inpos.preds90$f, positive = pos.preds90$f, date = seq(188,277,1))
  1. Predictions
#7 day predictions
#us7 = predict(us.mlr.modelfit, newxreg = regressions7)
#90 day predictions
#us90 = predict(us.mlr.modelfit, newxreg = regressions90)
  1. Plotted Forecasts
#plot(seq(1,187,1), us$hospitalizedCurrently, type = "l", xlim = c(0,195), ylim = c(0,475000), ylab = "Currently Hospitalized COVID Cases", main = "7 Day Forecast for COVID Hospitalized Cases")
#lines(seq(188,194,1), us7$pred, type = "l", col = "red")
  1. ASE
#usSmall = us[1:150,]
#us.ksfit = lm(hospitalizedCurrently~positive+)


#EXAMPLE CODE
#Find ASE  Need to forecast last 30 of known series.  
#CMsmall = CM[1:478,]
#ksfit = lm(cmort~temp+part+Week, data = CMsmall)
#phi = aic.wge(ksfit$residuals)
#fit = arima(CMsmall$cmort,order = c(phi$p,0,0), seasonal = list(order = c(1,0,0), period = 52), xreg = cbind(CMsmall$temp, CMsmall$part, CMsmall$Week))

#AIC(fit) #2990.101

#last30 = data.frame(temp = CM$temp[479:508], part = CM$part[479:508], CMWeek = seq(479,508,1))
#get predictions
#predsCMort = predict(fit,newxreg = last30)

#ASE = mean((CM$cmort[479:508] - predsCMort$pred)^2)
#ASE #84.2398

#plot(seq(1,508,1), CM$cmort, type = "l",xlim = c(0,528), ylab = "Cardiac Mortality", main = "Last 30 Week Cardiac Mortality Forecast")
#lines(seq(479,508,1), predsCMort$pred, type = "l", col = "red")
  1. Rolling Windowed ASE

VAR Model

Since we are working with a seasonal and transformation component but it requires an additional transformation for the total positive COVID cases to become stationary for evaluation, we will only use positive increase of cases to predict currently hospitzliaed cases for the VAR model.

  1. Differenced and Seasonal Transformation VAR Model
#Positive Cases Transformations
#pos.diff = artrans.wge(us$positive, phi.tr = 1, lag.max = 100)
#pos.diff.2 = artrans.wge(pos.diff, phi.tr = 1, lag.max = 100)
#pos.trans = artrans.wge(pos.diff.2,phi.tr = c(0,0,0,0,0,0,1), lag.max = 100)
#Increase Positive Transformation
inpos.diff = artrans.wge(us$positiveIncrease, phi.tr = 1, lag.max = 100)

inpos.trans = artrans.wge(inpos.diff,phi.tr = c(0,0,0,0,0,0,1), lag.max = 100)

#Current Hospitalization Transformations
us.diff = artrans.wge(us$hospitalizedCurrently, phi.tr = 1, lag.max = 100)

currhosp.trans = artrans.wge(us.diff,phi.tr = c(0,0,0,0,0,0,1), lag.max = 100)

  1. Use Var select to find best model and fit model: p = 8 for lowest AIC
#VARselect to choose lag
VARselect(cbind(currhosp.trans, inpos.trans), lag.max = 10, type= "both")
## $selection
## AIC(n)  HQ(n)  SC(n) FPE(n) 
##      8      7      7      7 
## 
## $criteria
##                   1            2            3            4            5
## AIC(n) 3.024021e+01 3.010098e+01 3.005507e+01 3.002311e+01 2.998028e+01
## HQ(n)  3.030034e+01 3.019117e+01 3.017532e+01 3.017342e+01 3.016066e+01
## SC(n)  3.038837e+01 3.032322e+01 3.035139e+01 3.039351e+01 3.042476e+01
## FPE(n) 1.358825e+13 1.182271e+13 1.129309e+13 1.093932e+13 1.048282e+13
##                   6            7            8            9           10
## AIC(n) 2.989906e+01 2.969862e+01 2.969858e+01 2.973472e+01 2.971028e+01
## HQ(n)  3.010951e+01 2.993912e+01 2.996915e+01 3.003536e+01 3.004098e+01
## SC(n)  3.041763e+01 3.029126e+01 3.036530e+01 3.047553e+01 3.052516e+01
## FPE(n) 9.667838e+12 7.914804e+12 7.918312e+12 8.214745e+12 8.022340e+12
#fit model
usdiff.var.fit = VAR(cbind(currhosp.trans, inpos.trans), type = "both", p = 2)
  1. Predictions for Difference
#7 day predictions
pred.var7 = predict(usdiff.var.fit, n.ahead = 7)
pred.var7
## $currhosp.trans
##           fcst     lower    upper       CI
## [1,] -119.8417 -2538.900 2299.217 2419.059
## [2,] -439.6708 -2909.280 2029.938 2469.609
## [3,] -171.5330 -2832.443 2489.377 2660.910
## [4,] -215.1457 -2895.958 2465.666 2680.812
## [5,] -146.9248 -2859.669 2565.820 2712.745
## [6,] -147.8477 -2867.280 2571.585 2719.433
## [7,] -120.3876 -2846.115 2605.340 2725.727
## 
## $inpos.trans
##              fcst     lower    upper       CI
## [1,]  795.5685172 -4387.729 5978.866 5183.297
## [2,]    0.8775186 -5671.247 5673.002 5672.125
## [3,] -389.3039566 -6073.653 5295.045 5684.349
## [4,] -116.4012417 -5816.574 5583.771 5700.172
## [5,] -117.3898871 -5820.079 5585.299 5702.689
## [6,] -184.3734040 -5887.147 5518.400 5702.774
## [7,] -146.9481382 -5850.038 5556.142 5703.090
#90 day predictions
pred.var90 = predict(usdiff.var.fit, n.ahead = 90)
pred.var90
## $currhosp.trans
##             fcst     lower    upper       CI
##  [1,] -119.84165 -2538.900 2299.217 2419.059
##  [2,] -439.67081 -2909.280 2029.938 2469.609
##  [3,] -171.53298 -2832.443 2489.377 2660.910
##  [4,] -215.14573 -2895.958 2465.666 2680.812
##  [5,] -146.92480 -2859.669 2565.820 2712.745
##  [6,] -147.84770 -2867.280 2571.585 2719.433
##  [7,] -120.38757 -2846.115 2605.340 2725.727
##  [8,] -118.37251 -2846.017 2609.272 2727.645
##  [9,] -108.12040 -2837.120 2620.879 2729.000
## [10,] -105.73962 -2835.244 2623.765 2729.504
## [11,] -102.15977 -2831.968 2627.649 2729.808
## [12,] -101.08856 -2831.024 2628.847 2729.936
## [13,]  -99.99271 -2829.999 2630.013 2730.006
## [14,]  -99.85993 -2829.897 2630.178 2730.037
## [15,]  -99.87055 -2829.924 2630.183 2730.054
## [16,] -100.25788 -2830.319 2629.804 2730.062
## [17,] -100.76737 -2830.833 2629.298 2730.065
## [18,] -101.43558 -2831.503 2628.632 2730.067
## [19,] -102.17597 -2832.244 2627.892 2730.068
## [20,] -102.98767 -2833.056 2627.081 2730.069
## [21,] -103.83823 -2833.907 2626.231 2730.069
## [22,] -104.72174 -2834.791 2625.347 2730.069
## [23,] -105.62530 -2835.694 2624.444 2730.069
## [24,] -106.54446 -2836.614 2623.525 2730.069
## [25,] -107.47370 -2837.543 2622.595 2730.069
## [26,] -108.41045 -2838.480 2621.659 2730.069
## [27,] -109.35220 -2839.421 2620.717 2730.069
## [28,] -110.29758 -2840.367 2619.772 2730.069
## [29,] -111.24544 -2841.315 2618.824 2730.069
## [30,] -112.19506 -2842.264 2617.874 2730.069
## [31,] -113.14590 -2843.215 2616.923 2730.069
## [32,] -114.09760 -2844.167 2615.972 2730.069
## [33,] -115.04990 -2845.119 2615.019 2730.069
## [34,] -116.00261 -2846.072 2614.067 2730.069
## [35,] -116.95563 -2847.025 2613.114 2730.069
## [36,] -117.90884 -2847.978 2612.160 2730.069
## [37,] -118.86220 -2848.931 2611.207 2730.069
## [38,] -119.81566 -2849.885 2610.253 2730.069
## [39,] -120.76919 -2850.838 2609.300 2730.069
## [40,] -121.72277 -2851.792 2608.346 2730.069
## [41,] -122.67638 -2852.746 2607.393 2730.069
## [42,] -123.63002 -2853.699 2606.439 2730.069
## [43,] -124.58368 -2854.653 2605.485 2730.069
## [44,] -125.53734 -2855.606 2604.532 2730.069
## [45,] -126.49102 -2856.560 2603.578 2730.069
## [46,] -127.44470 -2857.514 2602.624 2730.069
## [47,] -128.39838 -2858.468 2601.671 2730.069
## [48,] -129.35207 -2859.421 2600.717 2730.069
## [49,] -130.30576 -2860.375 2599.763 2730.069
## [50,] -131.25945 -2861.329 2598.810 2730.069
## [51,] -132.21314 -2862.282 2597.856 2730.069
## [52,] -133.16683 -2863.236 2596.902 2730.069
## [53,] -134.12053 -2864.190 2595.949 2730.069
## [54,] -135.07422 -2865.143 2594.995 2730.069
## [55,] -136.02791 -2866.097 2594.041 2730.069
## [56,] -136.98161 -2867.051 2593.088 2730.069
## [57,] -137.93530 -2868.004 2592.134 2730.069
## [58,] -138.88899 -2868.958 2591.180 2730.069
## [59,] -139.84269 -2869.912 2590.226 2730.069
## [60,] -140.79638 -2870.866 2589.273 2730.069
## [61,] -141.75007 -2871.819 2588.319 2730.069
## [62,] -142.70377 -2872.773 2587.365 2730.069
## [63,] -143.65746 -2873.727 2586.412 2730.069
## [64,] -144.61116 -2874.680 2585.458 2730.069
## [65,] -145.56485 -2875.634 2584.504 2730.069
## [66,] -146.51854 -2876.588 2583.551 2730.069
## [67,] -147.47224 -2877.541 2582.597 2730.069
## [68,] -148.42593 -2878.495 2581.643 2730.069
## [69,] -149.37963 -2879.449 2580.690 2730.069
## [70,] -150.33332 -2880.402 2579.736 2730.069
## [71,] -151.28701 -2881.356 2578.782 2730.069
## [72,] -152.24071 -2882.310 2577.828 2730.069
## [73,] -153.19440 -2883.264 2576.875 2730.069
## [74,] -154.14809 -2884.217 2575.921 2730.069
## [75,] -155.10179 -2885.171 2574.967 2730.069
## [76,] -156.05548 -2886.125 2574.014 2730.069
## [77,] -157.00918 -2887.078 2573.060 2730.069
## [78,] -157.96287 -2888.032 2572.106 2730.069
## [79,] -158.91656 -2888.986 2571.153 2730.069
## [80,] -159.87026 -2889.939 2570.199 2730.069
## [81,] -160.82395 -2890.893 2569.245 2730.069
## [82,] -161.77765 -2891.847 2568.291 2730.069
## [83,] -162.73134 -2892.800 2567.338 2730.069
## [84,] -163.68503 -2893.754 2566.384 2730.069
## [85,] -164.63873 -2894.708 2565.430 2730.069
## [86,] -165.59242 -2895.662 2564.477 2730.069
## [87,] -166.54612 -2896.615 2563.523 2730.069
## [88,] -167.49981 -2897.569 2562.569 2730.069
## [89,] -168.45350 -2898.523 2561.616 2730.069
## [90,] -169.40720 -2899.476 2560.662 2730.069
## 
## $inpos.trans
##               fcst     lower    upper       CI
##  [1,]  795.5685172 -4387.729 5978.866 5183.297
##  [2,]    0.8775186 -5671.247 5673.002 5672.125
##  [3,] -389.3039566 -6073.653 5295.045 5684.349
##  [4,] -116.4012417 -5816.574 5583.771 5700.172
##  [5,] -117.3898871 -5820.079 5585.299 5702.689
##  [6,] -184.3734040 -5887.147 5518.400 5702.774
##  [7,] -146.9481382 -5850.038 5556.142 5703.090
##  [8,] -153.2762415 -5856.392 5549.839 5703.115
##  [9,] -156.4642311 -5859.597 5546.669 5703.133
## [10,] -156.1862520 -5859.322 5546.949 5703.136
## [11,] -156.5209044 -5859.659 5546.618 5703.138
## [12,] -158.5053711 -5861.644 5544.634 5703.139
## [13,] -159.3557037 -5862.495 5543.784 5703.140
## [14,] -160.7370817 -5863.877 5542.403 5703.140
## [15,] -162.0720637 -5865.212 5541.068 5703.140
## [16,] -163.4544455 -5866.594 5539.686 5703.140
## [17,] -164.8303970 -5867.970 5538.310 5703.140
## [18,] -166.2496747 -5869.390 5536.890 5703.140
## [19,] -167.6623941 -5870.802 5535.478 5703.140
## [20,] -169.0893949 -5872.229 5534.051 5703.140
## [21,] -170.5190194 -5873.659 5532.621 5703.140
## [22,] -171.9533981 -5875.093 5531.187 5703.140
## [23,] -173.3894535 -5876.530 5529.751 5703.140
## [24,] -174.8277888 -5877.968 5528.312 5703.140
## [25,] -176.2670602 -5879.407 5526.873 5703.140
## [26,] -177.7073173 -5880.847 5525.433 5703.140
## [27,] -179.1481119 -5882.288 5523.992 5703.140
## [28,] -180.5893566 -5883.729 5522.551 5703.140
## [29,] -182.0308754 -5885.171 5521.109 5703.140
## [30,] -183.4726095 -5886.613 5519.667 5703.140
## [31,] -184.9144813 -5888.055 5518.226 5703.140
## [32,] -186.3564560 -5889.497 5516.784 5703.140
## [33,] -187.7984995 -5890.939 5515.342 5703.140
## [34,] -189.2405928 -5892.381 5513.899 5703.140
## [35,] -190.6827200 -5893.823 5512.457 5703.140
## [36,] -192.1248714 -5895.265 5511.015 5703.140
## [37,] -193.5670394 -5896.707 5509.573 5703.140
## [38,] -195.0092193 -5898.149 5508.131 5703.140
## [39,] -196.4514074 -5899.591 5506.689 5703.140
## [40,] -197.8936013 -5901.034 5505.246 5703.140
## [41,] -199.3357991 -5902.476 5503.804 5703.140
## [42,] -200.7779998 -5903.918 5502.362 5703.140
## [43,] -202.2202025 -5905.360 5500.920 5703.140
## [44,] -203.6624065 -5906.802 5499.478 5703.140
## [45,] -205.1046115 -5908.245 5498.035 5703.140
## [46,] -206.5468172 -5909.687 5496.593 5703.140
## [47,] -207.9890234 -5911.129 5495.151 5703.140
## [48,] -209.4312298 -5912.571 5493.709 5703.140
## [49,] -210.8734366 -5914.013 5492.267 5703.140
## [50,] -212.3156434 -5915.456 5490.824 5703.140
## [51,] -213.7578504 -5916.898 5489.382 5703.140
## [52,] -215.2000575 -5918.340 5487.940 5703.140
## [53,] -216.6422646 -5919.782 5486.498 5703.140
## [54,] -218.0844718 -5921.225 5485.056 5703.140
## [55,] -219.5266789 -5922.667 5483.613 5703.140
## [56,] -220.9688861 -5924.109 5482.171 5703.140
## [57,] -222.4110934 -5925.551 5480.729 5703.140
## [58,] -223.8533006 -5926.993 5479.287 5703.140
## [59,] -225.2955078 -5928.436 5477.845 5703.140
## [60,] -226.7377151 -5929.878 5476.402 5703.140
## [61,] -228.1799223 -5931.320 5474.960 5703.140
## [62,] -229.6221295 -5932.762 5473.518 5703.140
## [63,] -231.0643368 -5934.204 5472.076 5703.140
## [64,] -232.5065440 -5935.647 5470.634 5703.140
## [65,] -233.9487513 -5937.089 5469.191 5703.140
## [66,] -235.3909585 -5938.531 5467.749 5703.140
## [67,] -236.8331658 -5939.973 5466.307 5703.140
## [68,] -238.2753730 -5941.415 5464.865 5703.140
## [69,] -239.7175803 -5942.858 5463.422 5703.140
## [70,] -241.1597875 -5944.300 5461.980 5703.140
## [71,] -242.6019948 -5945.742 5460.538 5703.140
## [72,] -244.0442020 -5947.184 5459.096 5703.140
## [73,] -245.4864093 -5948.626 5457.654 5703.140
## [74,] -246.9286165 -5950.069 5456.211 5703.140
## [75,] -248.3708238 -5951.511 5454.769 5703.140
## [76,] -249.8130310 -5952.953 5453.327 5703.140
## [77,] -251.2552383 -5954.395 5451.885 5703.140
## [78,] -252.6974455 -5955.838 5450.443 5703.140
## [79,] -254.1396528 -5957.280 5449.000 5703.140
## [80,] -255.5818600 -5958.722 5447.558 5703.140
## [81,] -257.0240673 -5960.164 5446.116 5703.140
## [82,] -258.4662745 -5961.606 5444.674 5703.140
## [83,] -259.9084818 -5963.049 5443.232 5703.140
## [84,] -261.3506890 -5964.491 5441.789 5703.140
## [85,] -262.7928963 -5965.933 5440.347 5703.140
## [86,] -264.2351035 -5967.375 5438.905 5703.140
## [87,] -265.6773108 -5968.817 5437.463 5703.140
## [88,] -267.1195180 -5970.260 5436.021 5703.140
## [89,] -268.5617253 -5971.702 5434.578 5703.140
## [90,] -270.0039325 -5973.144 5433.136 5703.140
  1. Calculate Actual Forecasts from Predicted Differences
startingPoints7 = us$hospitalizedCurrently[181:187]
currHospForecasts7 = pred.var7$fcst$currhosp.trans[,1:3] + startingPoints7
startingPoints90 = us$hospitalizedCurrently[98:187]
currHospForecasts90 = pred.var90$fcst$currhosp.trans[,1:3] +startingPoints90
  1. Plotting Forecasts
#7 day Forecasts
plot(seq(1,187,1), us$hospitalizedCurrently, type = "l", xlim = c(0,195), ylab = "Currently Hospitalized COVID Patients", main = "7 Day Currently Hospitalized Patients Forecast")
lines(seq(188,194,1), as.data.frame(currHospForecasts7)$fcst, type = "l", col = "red")

#90 day Forecasts
plot(seq(1,187,1), us$hospitalizedCurrently, type = "l", xlim = c(0,280), ylab = "Currently Hospitalized COVID Patients", main = "90 Day Currently Hospitalized Patients Forecast")
lines(seq(188,277,1), as.data.frame(currHospForecasts90)$fcst, type = "l", col = "red")

  1. Calculate ASE #WHERE YOU LEFT OFF
#varASE = mean((us$hospitalizedCurrently[123:132]-pred.var$fcst$y1[1:10])^2)
#varASE